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Abstract 

In this paper we have explored computationally the solidification process of large nickel clus- 

ters. This process has the characteristic features of the first order phase transition occurring in a 

finite system. The focus of our research is placed on the elucidation of correlated dynamics of a 

large ensemble of particles in the course of the nanoscale liquid-solid phase transition through the 

computation and analysis of the results of molecular dynamics (MD) simulations with the corre- 
ct 

sponding theoretical model. This problem is of significant interest and importance, because the 
controlled dynamics of systems on the nanoscale is one of the central topics in the development of 
modern nanotechnologies. 

MD simulations in large molecular systems are rather computer power demanding. Therefore, 
in order to advance with MD simulations we have used modern computational methods based on 
the graphics processing units (GPU). The advantages of the use of GPUs for MD simulations in 
comparison with the CPUs are demonstrated and benchmarked. The reported speedup reaches 
factors greater than 400. This work opens a path towards exploration with the use of MD of a larger 
number of scientific problems inaccessible earlier with the CPU based computational technology. 
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I. INTRODUCTION 



The melting process in a macroscopically large system occurs at a certain temperature 
under fixed external pressure. This process is a first order phase transition and it manifests 
itself as a spike in the heat capacity of the system at the transition temperature. The 
reverse process of freezing or solidification occurs at a temperature below the thermodynamic 
melting point due to nucleation nature of the solidification phase transitions 1 . 

In small systems the finite volume that is available for the precursor formation alters the 
kinetics of the solidification phase transition. The study of supercooled metal droplets of 
the size of several micrometers was performed in early 1950s 2,3 . In Ref.- it was shown that 
the mercury droplets of 2 - 8 fim in diameter solidify at rates that are proportional to the 
droplet volume. In the systems of the nanometers size the ratio of the surface to volume 
atoms increases and the reduced binding energy of the surface atoms changes significantly 
the total energy of the system. The first thermodynamic model bridging the melting point of 
the clusters with their size was proposed by Pawlow more than a century ago^. In multiple 
subsequent works it was confirmed that the melting temperature of a small spherical particle 
decreases with the reduction of its radius 5- -. In small particles the relative fraction of the 
surface atoms is higher, which leads to the decrease of the melting temperature. This size 
effect have been confirmed for the clusters having diameters down to 2 nm 5 . Its explanation 
is based on the characteristic radial dependence of the ratio of the surface to volume energy 
of a finite system. 

However, for clusters having sizes smaller than 1-2 nm, the melting temperature is 
no longer a monotonic function of the cluster size. Experiments on sodium clusters Na^ 
with number of atoms N=50 - 360 have demonstrated that the melting temperature as a 
function of size shows a prominent irregular structure with local maxima^— . The origin of 
the nonmonotonic variation in the melting temperature with respect to cluster size lies in 
the interplay between electronic and geometric shell effects in the sodium clusters and the 
entropy change in the course of such nanoscale phase transitions 12 . 

Intensive theoretical efforts have been undertaken to identify the details of the geometric 
and electronic structures underlying the variations in the melting temperature^— . Experi- 
ments on small ion clusters of tin^i and gallium 2 ^ have confirmed the violation of the linear 
relationship between the reduction in the melting temperature and the inverse radius of the 
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cluster. It was discovered that the melting temperature of selected Sn^ and Ga^v clusters, 
of sizes N < 40, can considerably exceed the melting temperature of the corresponding bulk 
material" 1 ' 33 . This behavior was explained by the structural differences between the small 
clusters and the bulk^i. 

There are also other factors which affect the melting temperatures of small clusters. For 
instance, impurities play a significant role. Thus, in Ref.— it was demonstrated that a single 
atom impurity leads to significant changes in the thermodynamic properties of the Ni 147 
cluster. Also, the melting temperatures differ significantly in nanoalloys as compared to 
pure materials. In KefM it was shown that alloying iron clusters consisting of up to 2400 
atoms with carbon reduces their melting temperature by 100 - 150 K at a carbon concen- 
tration of 10% - 12%. Recently, the thermal behaviour of free and alumina-supported Fe-C 
nanoparticles has been investigated 25 . It was observed that the presence of the substrate 
raises the melting temperature of medium and large Fe(i_ x )ArC X Ar nanoparticles (x=0 - 0.16, 
N=80 — 1000) by 40 - 60 K— . In Ref.— the freezing-melting hysteresis associated with a 
free energy barrier between solid and liquid phases was investigated theoretically for the 
materials in pores. 

There are many papers devoted to the computer simulations of the melting process of 
metal clusters consisting of up to few hundreds of atoms. For the small systems one can 
utilize ab initio approaches based on Density Functional Theory 2 ^ or more empirical Car- 
Parinello 2 ^ or Tight-Binding schemes^. The aforementioned methods are capable of repro- 
ducing relatively accurately the structural and energetic properties of the clusters. However, 
computational complexity of ab initio methods does not allow for modelling of metal clusters 
consisting of thousands of atoms on the time scales exceeding nanoseconds. 

Investigations of melting of clusters of several thousands of atoms are carried out using 
empirical classical potentials for the description of interatomic interactions. The delocal- 
ized nature of <i-electrons in the transition metal clusters implies utilization of so-called 
many-body potentials, which account for two- and three-body interatomic interactions. The 
most widely used potentials for the description of interactions between nickels atoms are 
Finnis-Sinclair and Sutton-Chen potentials. In Ref.^2, MD simulations of melting and solidi- 
fication of Ni nanoclusters with up to 8007 atoms are performed with the use of the classical 
quantum-corrected Sutton-Chen potential^, and the characteristics of the liquid-solid phase 
transition are analysed for different cluster sizes. The investigation of melting and solidifi- 
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cation of Ni clusters consisting of up to 32000 atoms using the conventional Sutton-Chen 
potential is presented ivM. In that paper the reduction of the melting and solidification tem- 
peratures for the face-centred cubic metals is reported to be negatively correlated with the 
particle radius, and the Gibbs-Thomson coefficient is found to be proportional to the melt- 
ing point. In Ref.— the investigation of the effect of the cooling rate on the final structure 
of M06750 clusters is performed with the use of the Finnis-Sinclair potential. 

In the present paper we conduct MD simulations of Ni 20 57 clusters on the timescale up 
to 65 ns and discuss the kinetics of the solidification phase transition for the clusters as a 
function of the overcooling rate. From the analysis of MD simulations we extrapolate the 
thermodynamic properties for the clusters of arbitrary size up to the bulk. 

Focus of our research is placed on the elucidation of correlated dynamics of a large 
ensemble of particles in the course of the nanoscale liquid-solid phase transition. Obviously, 
this problem is of significant interest and importance, because the controlled dynamics 
of systems on the nanoscale is one of the central topics in the development of modern 
nanotechnologies. For the purposes of this analysis we have chosen the melting process of 
large nickel clusters and performed a systematic theoretical analysis of their dynamics in 
the course of melting and solidification. 

The choice of Ni clusters for these studies is motivated by their high chemical and catalytic 
reactivity, unique properties, and multiple applications in nanostructured materials^. An 
important example of such an application is the process of the catalytically activated growth 
of carbon nanotubes. The thermodynamic state of the catalytic nanoparticle plays a crucial 
role in the carbon nanotube growth^. The important question is whether the catalytic 
nanoparticle is molten or frozen during the nanotube growth process. It was demonstrated 
that when carbon nanotubes are grown on large 3-4 nm iron nanoparticles at temperatures 
lower than 1200 K, the catalytic particle is not completely molten™. Thus, the mechanism 
of nanotube growth can be governed by the surface melting of the cluster. Therefore, the 
advanced MD simulations of melting of large Ni clusters are important for a reliable eval- 
uation of the conditions at which the carbon nanotube growth process takes place^ 1 ^ 1 ^ 1 ^. 
In order to achieve this task it is necessary to perform MD simulations for relatively large 
nanoparticle sizes and large simulation times, which imposes the use of the most advanced 
computational techniques based on the GPUs. The advantages of this technology for MD 
simulations in comparison with the CPU based one is demonstrated and benchmarked. The 
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reported speedup reaches factors greater than 400. This work opens a path towards ex- 
ploration of a larger number of scientific problems inaccessible earlier with the CPU based 
computational technology. 

The CPU calculations of the nickel clusters were carried out using a multi-purpose com- 
puter code MesoBioNano Explorer (MBN Explorer)^. MBN Explorer allows to use a broad 
variety of interatomic potentials, to model different molecular systems, such as atomic clus- 
ters, fullerenes, nanotubes, polypeptides, proteins, DNA, composite systems, nanofractals, 
etc. Despite the universality, the computational efficiency of MBN Explorer is comparable 
(and in some cases even higher) than the computational efficiency of other software packages. 

For the purposes of this work we have adopted some of the MD algorithms of MBN 
Explorer to run them on GPUs. In particular, we have rewritten the part of the code 
responsible for calculations of the Sutton-Chen potential and force using Open Computing 
Language (OpenCL). OpenCL is a framework for writing programs that execute across 
heterogeneous platforms consisting of CPUs, GPUs, and other processors. The details of 
the implementation are discussed in the following section. 

The paper is organized as follows. In section [Til we introduce the Sutton-Chen potential, 
describe the theoretical model of nucleation and growth of solid state precursors in the course 
of the solidification phase transition, and present the details of the computational approach 
utilized in the work. In section IHII we demonstrate the results of MD simulations of Ni2047 
clusters and analyse radial distribution function, diffusion coefficients for molten and solid 
states of the cluster, and the correspondence of MD simulation results with the theoretical 
model for the solidification rate. In section ITVl we draw conclusions to the paper. 

II. THEORETICAL AND COMPUTATIONAL METHODS 
A. Interaction potential for nickel atoms 

The study of structural and dynamical properties of transition-metal clusters is a chal- 
lenging task due to the presence of unfilled valence d orbitals. The high density of the 
d states and their delocalized character make the direct ab initio methods computation- 
ally very demanding for clusters larger than several tens of atoms^. In order to describe 
the structure of clusters of larger sizes, one needs to use approximate methods and model 
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interatomic potentials. 

An effective approach for the study of transition-metal clusters is the embedded-atom 
method 4 ^"— , which takes into account many-body effects. The latter appears through the 
inhomogeneous electron density of the system. In this paper, an MD study of nickel clusters 
has been performed using the Sutton-Chen 4 ^ many-body potential, which belongs to the 
family of the embedded-atom types of potentials. The Sutton-Chen potential 4 ^ has been 
shown to reproduce bulk and surface properties of transition metals and their alloys with 
sufficient accuracy see, e.g. Refs.— £t£L and references therein. The applicability of the 
Sutton-Chen potential 44 to Ni clusters has been proven by the direct comparison of the 
optimized structures and the binding energies of small Ni clusters obtained within ab initio 
method and with the use of the Sutton-Chen potential^. 

The potential energy of the finite system within the Sutton-Chen model has the following 
form: 



Here is the distance between atoms i and j, e is a parameter with dimension of energy, 
a is the lattice constant, c is a dimensionless parameter, and n and m are positive integers 
with n > m. The parameters provided by Sutton and Chen for nickel have the following 
values:^ e = 1.5707 • 1(T 2 eV, a = 3.52 A c = 39.432, n = 9, and m = 6. 



B. Theoretical model of the solidification kinetics 

In the liquid-solid phase transition the process of the formation of the new crystal phase 
is initiated by the formation of the stable precursor of the solid state. The energy of the 
formation of the precursor E p (r) of radius r can be written as follows^: 




(1) 



where 




(2) 
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E P (r) = — (3a - rAf) , 



(3) 



where a is the liquid-solid surface tension coefficient and Af is the difference of the free 
energy densities for the liquid and the solid states. Here we consider the process of crystal- 
lization from the overcooled liquid phase, therefore Af is positive and can be evaluated as 
follows: 



Af = fi - f s , (4) 

where fi and f s are the free energy densities for the liquid and solid phases, correspondingly. 
fi and f s can be written as: 

3NkT 

fi = e i0 + — Ts u (5) 

Vcl 

f i 3NkT T 

fs = e s0 + — Ts s , (6) 

Vcl 

where e/ , e s0 , sj, s s , N and V c \ are the ground-state energy densities of the liquid and solid 
phases, entropies densities of the liquid and solid states, number of atoms in the cluster 
and cluster volume, correspondingly. The factor 3kT accounts for the energy stored in the 
thermal vibrations of the atoms at finite temperature. Using Eqs. (JMH]), Eq. (j4]) can be 
written as follows: 

Af = As(T-T ), (7) 

where As = s s — si and T is the phase transition temperature. For the derivation of Eq. ([7]) 
we have accounted for the fact that free energy densities for both phases are equal at the 
phase transition temperature, i.e. e s o — T s s = eio — T si. 



Eq. ([3]) has a maximum at certain value of r, corresponding to the critical size r c of the 
precursor of the solid phase. If the size of the precursor is larger than r c its further growth 
will reduce the free energy of the system, and the crystal phase will expand. The value of 
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r c corresponds to the size at which the derivative over r of the Eq. ([3]) is equal to zero, i.e. 

(8) 



2a 



Substituting r c from Eq. (jSJ) into Eq. (j3J) one obtains the following value of the energy 
associated with the formation of the precursor of critical size: 

167TOJ 3 



E p {r c ) 



3As 2 (T - T 



,2 ' 



(9) 



where we have used Eq. ([7]) to substitute the expression for Af. 

Note that the precursors of the crystal phase are very improbable to be formed in the 
vicinity of the cluster surface, since the surface atoms are in a liquid-like state even at the 
temperatures when the core of the cluster is in the solid state. 

The rate of the formation of the precursors of the solid phase can be calculated as follows: 



k = A exp 



lQna 3 



3k B TAs 2 (T — T ) j 

where k B is the Boltzmann constant and A is the precursor formation frequency. 



(10) 



Computational approach 

The numerical solution of the MD equations becomes a challenging task if one studies 
systems with a large number of atoms. Generally, the crucial time-consuming part of the 
computation is the calculation of the forces between the particles at each time step. For 
two-body forces this implies a computational effort rising with N 2 , where N is the number 
of atoms, see Fig. [TJ An implicit many-body structure of the force as it is generated by the 
second term in the Sutton-Chen potential, Eq. (JT]), can potentially increase the calculation 
effort. 



As the force calculation essentially implies the same type of calculation repeated many 
times over, this task is ideally suited to be implemented using a parallel-programming ap- 
proach. If there are only few different types of forces involved, as is the case discussed 
here, a single-instruction multiple-data (SIMD) hardware environment can be adopted in a 
natural way, since the same instructions are used for all pairs of particles. Graphics cards 
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Number of particles (in thousands) 

Figure 1. Dependence of computation time on the system size. From Fig. it is seen that the 
time is proportional to the squared number of particles, i.e. 0(N 2 ). Note that the time per step 
is multiplied by a thousand for the GPU case. 

(GPUs) fulfill these requirements in an ideal way. We therefore ported our original CPU 
based code MBN Explorer 38 for CUDA as well as OpenCL programming frameworks. 

We reduced the many-body part of the calculation to two successive two-body problems, 
in the first sweep calculating all pi in Eq. (J2J) and then determining the forces on the atoms. 

The calculations presented here were performed on the LOEWE compute cluster at Frank- 
furt University, which features 778 graphics cards Radeon HD 5870. Each card includes 1600 
streaming processors clocked at 850 MHz. The results GPU calculations were compared to 
a single core of Intel Core i7 870 CPU. 

Fig.|2]shows the speed-up of a MD simulation of nickel clusters of different sizes comparing 
the CPU and GPU versions of the code. As can be seen one can achieve speed-ups of more 
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Number of particles (in thousands) 

Figure 2. Speed-up of the computational time for GPU as compared to the CPU version. 

than two orders of magnitude for larger clusters. These results are consistent with theoretical 
estimation of performance of GPU and CPU in single precision. AMD Radeon specified to 
have 2720 GFLOPS in single precision and one core of Core i7 CPU can be estimated as 
nearly 15 GFLOPS in double precision. In double precision the difference would be less as 
AMD Radeon specified to have 544 GFLOPS in double precession. It means that for double 
precision code one can expect speed-up of nearly one hundred over single threaded code. 

In our code each particle was assigned to a certain thread. It implies that if the number of 
particles is smaller then the number of stream processors GPU will not be completely loaded. 
In order to utilize the GPU efficiently the number of particles should be an integer multiple 
of the number of the stream processors. This can be seen on speed-up plots presented in 
Figs. [U EJ Speed-up growth is nearly linear until 10 thousands of particles. For the number 
of particles exceeding 10 thousands the speedup saturates to a constant level, corresponding 
to full utilization of the GPU. In the presented benchmarks the GPU computation time was 
averaged over one thousand iterations in order avoid accounting for random slowdowns and 
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lags. CPU computation time was averaged over 10 iterations since random slowdowns of 
CPU are expected to have minor influence on the computational time. 

Current calculations were performed as full N-body calculation without spatial cut-off 
for the interaction of the atoms, i.e., taking into account interactions between all atoms of 
the cluster. Therefore, in this work we have studied the cluster consisting of ~2000 atoms, 
where a cut-off does not have a substantial advantage in computing time, but the speedup 
due to utilization of GPUs is about 100. The efficient implementation of interaction cutoffs 
in the GPU code is in progress and will be adopted and discussed in forthcoming work. 

The MD simulations were carried out using a Verlet integrator with a time-step of 1 fs. 
The temperature control was carried out using a Langevin thermostat with a damping 
constant of 10 ps _1 . 

III. RESULTS AND DISCUSSION 
MD simulations of Ni2047 clusters 

We have performed MD simulations of M2047 clusters. The initial geometry of the cluster 
was chosen to resemble icosahedral symmetry. Then the particle was exposed to the heating- 
cooling cycle with the heating/cooling rate of 1 K/ps in the temperature range from 400 to 
1400 K. From Fig. [3] is is seen that the final structure of the crystallized particle (blue dots 
at low temperatures) has lower energy than the initial icosahedral structure. The method 
to obtain molecular structures with higher binding energy by exposing the heated system 
to the cooling temperature bath is widely used technique in MD simulations and is known 
as simulated annealing^. The lowering of the nickel cluster's energy after resolidification 
form icosahedral structure happens due to the formation of regions with FCC symmetry in 
the relatively large nickel clusters^. From Fig. [3] it is also seen a prominent hysteresis in 
the melting and crystallisation phase transitions due to the finite speed of the heating and 
cooling rates. Note that the width of the hysteresis is almost 400 K for the heating/ cooling 
rate of 1 K/ps, which shows that the crystallization process in the system is associated with 
transfer over a relatively high free energy barrier. 

In the following section we investigate the dynamics of the crystallisation transition. 
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Temperature K 

Figure 3. Dependence of the total energy per atom in heating (red) and cooling (blue) simulations. 
In the insets are shown the structures of N12047 in the molten (top-left) and solidified (bottom-right) 
states. The fragment of crystal structure is shown by blue spheres inside the crystallized structure. 

A. Solidification process in Ni2047 clusters 

By means of MD simulations we have investigated the solidification kinetics of Ni 20 47 
cluster as a function of the amount of overcooling. The initial structure of the cluster was 
taken from MD simulations above the melting temperature. Then the initially molten cluster 
was simulated at various temperatures below the melting point. We have investigated the 
time needed for the system to change its phase from the molten to the solid state as a function 
of the temperature below the phase transition point. The simulations were performed at 
740 K, 760 K, 780 K, 800 K, 820 K, 830 K and 835 K. For each temperature up to 15 
independent simulations were produced. Fig. 0] shows the dependence of the total energy 
of the system on the simulation time for 5 randomly chosen trajectories for each value of 
temperature in the range between 740 and 820 K. Note the logarithm scale on the horizontal 



12 



axes. 



-3,88- 



> 

0) 



! -740K 
T=800-K 



111*1 i < i i i * i 1 1 

- T-760K ■ TWBDK 
T=920K 




I -4,00 

to 

o 



-4,04 



10' 



* i i i i 1 1 1 i 

10 3 
Time (ps) 




i p i i i i | 



10' 



Figure 4. Dependencies on time of the total energy per atom for the cluster being initially in the 
molten state at different temperatures of the thermostat. Solidification phase transition occurs 
faster in the systems with lower temperature. 



From Fig. H] it is seen that at certain moments of time the energy of the system abruptly 
changes. These moments correspond to the transition of the system from the molten to the 
solid state. Indeed, the solid state of the cluster corresponds to the increased binding energy 
between the atoms due to the formation of the regular crystal lattice, which leads to the 
lowering of the total energy of the system. 



B. Radial distribution function as an indicator of the phase transition 

In this section we analyse to which extent the radial distribution function (RDF) for 
atoms in the cluster becomes affected by the solidification phase transition. RDF is defined 
as: 
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g (r) = p(r)/p, (11) 

where p(r) is the density of atoms at distance r form a reference particle and p is an average 
density. g(r) characterises the number of particles at the certain radial distance from a 
reference particle. 




Radius (A) 

Figure 5. Radial distribution function for the most central atoms of N 22047 cluster at different 
instances of simulation time. For the selected for analysis trajectory the solidification transition 
had occurred at t pt = 5 ns. 

We have calculated the RDF for nickel clusters in solid and molten states. The following 
procedure was adopted for the calculation of the RDF. At certain instance of time the atoms 
located within 3 A of the cluster center of mass were selected. For these atoms the the RDF 
was calculated with a distance bin of 0.4 A and averaged over 1 consequent nanosecond 
of simulation, which results in averaging over 100 cluster structures, since we have written 
the structure of the system each 10 ps of simulation. The RDF for the central atoms was 
calculated as an average over RDFs for atoms of each of 100 structures. 

We analysed one particular trajectory of the MD simulations at the thermostat tem- 
perature equal to 820 K (see orange trajectories in Fig. H]). For the chosen trajectory the 
solidification transition occurs 5 ns after the start of the simulation. The RDF for that 
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trajectory at various instances of time is plotted using B-splines in Fig. [5j 

From Fig. [5] it is seen that RDFs at 1,5 and 3 ns (black and red curves) are substantially 
different form those plotted for later time instances. This, obviously, can be explained by 
the formation of cluster crystalline structure upon solidification. Indeed, peaks at 5 and 
7 A are more pronounced for the RDFs at 6 .. 9 ns, since the crystalline structure has a 
long-range order. Note that RDF at 4,5 ns attain many features of RDFs characteristic 
for the crystallized state, however the crystallization transition takes place at 5 ns. This 
is related to the fact the calculated RDF was averaged over 1 ns in order to increase the 
statistics. Therefore, RDF at 4,5 ns was partially constructed from the cluster structures 
being in crystalline state. 

Despite that RDF can be used for the characterisation of the melting-solidification tran- 
sition in the system, it is not entirely clear how to define uniquely the phase transition 
moment from RDF analysis, especially when RDFs are calculated with a reduced amount 
of sampling data, which is often the case for finite systems. In the next subsection we re- 
port on the analysis of the diffusion coefficient behaviour which turns out to be much more 
convenient quantity for the characterization of the liquid-solid phase transition. 

C. Diffusion coefficient as a fingerprint of solidification phase transition 

The diffusion coefficient of a particle is defined as follows: 



D = £g, (12) 
2zAt ' v ; 

where (Ar 2 ) is a mean-square displacement of a particle per time At, and z is dimensionality 
of space (3 for 3-dimensional diffusion) 52 . 

We have calculated the self-diffusion coefficient of nickel atoms located in the central part 
of the cluster at different instances of time. For the analysis we have chosen 3 out five MD 
trajectories of the A^" 2 o57 cluster conducted with a thermostat temperature equal to 820 K. 
The total energy of the system and the diffusion coefficient as a functions of time are plotted 
for each trajectory in Fig. [6j 

The following procedure was used for the calculation of the self-diffusion coefficient: the 
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Figure 6. Left axis: dependencies of the total energy per atom on time. Colours indicate different 
MD trajectories. Right axis: dependence of diffusion coefficient on time. 

part of the MD trajectory starting 0,5 ns before and ending 0,5 ns after each given instant of 
time was selected and the coordinates of all atoms were recorded. For each structure cluster 
was translated and rotated as a whole in order to minimize the root-mean-square displace- 
ment of atoms from a reference structure. This procedure allows to exclude rotational and 
translational motion of the cluster when calculating the diffusion coefficient. Then the se- 
lected part of the trajectory was split into 10 segments and the diffusion coefficient was 
calculated for all atoms located within 10 A of the cluster center of mass using Eq. ffT2l . 
The final value of the diffusion coefficient was calculated as an average over all central atoms 
and 10 trajectory segments. Since the diffusion coefficient was calculated as an average over 
1 ns of the MD simulation, in Fig. [6] we show the corresponding horizontal error bars for 
each value of the diffusion coefficient. 

From Fig. |6] it is seen that for all 3 trajectories the diffusion coefficient drops abruptly 
at certain moments of time exactly corresponding to the moment of the solidification phase 
transition. The value of the diffusion coefficient for the molten state is ~ 2.3 x 10 -6 cm 2 /s. 
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The diffusion coefficient in the solid state appears to be smaller than 10 -5 .. 10 -6 cm 2 /s. We 
were not able to evaluate it accurately for a given set of data since it requires more than 5 
ns to observe even single event of exchange of mutual positions for a pair of atoms at 820 K. 

From Fig. [6] it is seen that the dependence of the diffusion coefficient on time qualita- 
tively reproduces well the time dependence of the total energy. One can use the diffusion 
coefficient as a reliable quantity to characterize the solidification phase transition in the 
system especially if calculation of the total energy of the system during MD simulations is 
undesirable, which might be the case for certain GPU-based computer codes. 



D. Kinetics of the solidification phase transition 

In this subsection we evaluate thermodynamic characteristics of the solidification phase 
transition on the basis of the analysis of kinetics of this process obtained from MD simula- 
tions. 

We calculate the latent heat of the phase transition as a difference of the averaged total 
energy of the system in molten and solid states. For the NI2047 cluster the latent heat of the 
phase transition AE is equal to 127.8 eV. The entropy density change between solid and 
molten states of the cluster can be calculated as follows: 

A, = ^, (13) 

where To and V c \ is a phase transition temperature and cluster volume, correspondingly. 
Note, that here we neglect the fact that the surface layer of the atoms of the cluster remains 
liquid after the solidification of the cluster core. 

In Table H] we present the life time of the system in overcooled state for temperature range 
between 780 K and 835 K. For our analysis we do not take into account simulations performed 
for temperatures lower than 780 K since the lifetime of the clusters in the overcooled state 
for that temperatures is comparable with the time necessary for the system to exchange its 
energy with the thermostat in the course of the solidification process (see Fig. H]). The time 
of the phase transition was defined as a moment at which the energy of the system starts to 
decrease rapidly due to the formation of the solid phase. The last row in Table [J shows the 
escape rate from the overcooled states, which is defined as: 
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Figure 7. Fraction of M2057 clusters in the molten state as a function of simulation time for 
various temperatures below the phase transition point. Dots show the results obtained from MD 
simulations. Numerical fit of the results using Eq. (I10p is shown by solid lines. 



K=(<t>y\ (14) 

where < t > denotes averaging of the lifetime over the simulation runs for the corresponding 
temperature. 

The fraction of clusters in the liquid state below the phase transition temperature as a 
function of time, c(t), can be written as follows: 

c(t) = exp (—Kt) , (15) 
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Figure 8. Solidification rate of M2057 cluster for various temperatures below the phase transition 
point. Error bars correspond to a single standard deviation calculated for a set of data for each 
temperature value. Numerical fit of the results using Eq. (fT6|) is shown by solid line. 

where k can be evaluated using Eq. (JHJ). In FigJB] the results of the MD simulations are 
shown by dots. Solid lines correspond to the exponential reduction of the fraction of liquid 
clusters as derived using Eq. (|T5|) . 

Knowing the dependence of the solidification rate of the clusters on temperature we can 
fit the results of MD simulations to the Eq. (ITUi) as follows: 



ln{K - lfs) = a 'k B T(T-T r (16) 

where a, b and T are fitting parameters. The result of the fit of the MD simulations is shown 
in Fig. [BJ and the fitting parameters have the following values: a=1.24, 6=31053 eV-K 2 
and To=1029 K. Therefore, according to the results of the fitting, the precursor formation 
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3320 


3 


230 


180 


2260 


9340 


5890 


4 


240 


550 


3600 


18130 


14680 


5 


280 


820 


5020 


19080 


20008 


6 


310 


890 


5130 


19760 


22240 


7 


360 


1260 


6670 


22640 


28640 


8 


370 


1320 


7430 


26670 


30230 


9 


380 


1390 


7800 


27670 


41670 


10 


530 


1600 


9450 


36380 


63180 


11 


540 


1830 


11620 


— 


— 


12 


570 


1930 


23400 






13 


740 


2210 








14 


750 


2800 








15 


2000 


3610 








mean (ps) 
re (10- 5 ps _1 ) 


503 
199 


1378 
72.6 


6440 
15.5 


18731 
5.34 


23102 
4.33 



Table I. Life time in picoseconds of the overcooled liquid state for different MD trajectories at 
various temperatures. 15 simulations were performed for 780 K and 800 K, 12 for 820 K and 10 
simulations for 830 K and 835 K. 



frequency (coefficient A in Eq. (llOjl ) is 3.5 ps -1 and phase transition temperature of the 
bulk nickel material for the used parametrization of the Sutton-Chen potential is equal to 
1029 K. This value of phase transition temperature is lower than that, reported in^, where 
the phase transition of the bulk nickel for the similar parametrization of the Sutton-Chen 
potential was calculated to be 1160 K. The discrepancies between our results and the results 
of previous calculation are largely because of limited statistics for the calculation of the 
overcooled clusters lifetimes. 

Evaluating As using Eq. ( f!3|) one can calculate the liquid-solid surface tension coefficient 
a as a real root of the following equation: 

Solving Eq. ( fTTI) one obtains a = 0.53 eV/nm 2 or 85 mJ/m 2 . This value of a is approximately 
two times lower than that reported in Ref.— for Nickel crystal-liquid tension coefficient. The 
discrepancy can be attributed to the various parametrization of the Sutton-Chen potential 
used in our and that works. 
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Knowing the surface tension coefficient a one can calculate the critical size of the precursor 
of the solid phase for different temperatures using Eq. (jSJ) and Eq. f|T3|) . From that equations 
it follows that the critical precursor consists of 26, 33, 41, 53, 70, 80 and 87 atoms for the 
temperature equal to 740, 760, 780, 800, 820, 830 and 835 K, correspondingly. 

From Eq. (jSj) it is seen that the critical size of the precursor grows when the temperature 
of the system approaches the phase transition temperature of the bulk. However, for finite 
systems the size of the critical precursor can not be larger than the cluster size. We can use 
this idea to calculate the phase transition temperature of the finite size cluster systems as 
follows: 



where N denotes the number of particles in the cluster and p is the particle density. This 
equation is also known as Gibbs-Thomson equation describing the suppression of the melting 
point of the solid particles in their own fluid. From Eq. ( TT8|) it follows that the suppression 
of the melting point is inversely proportional to the number of atoms in the cluster in the 
power of one third, and the melting point of the cluster consisting of 2057 atoms is equal to 
962 K. 

IV. CONCLUSIONS 

In this paper we have conducted classical MD simulations of the M2047 cluster with the 
use of many-body Sutton-Chen potential on the time scales up to 65 ns. For the purposes 
of this work we have developed an efficient software code capable of performing large-scale 
MD simulations on GPUs. We have demonstrated that with the use of GPU technology 
one can achieve the speed-up of computations up to 400 times as compared to single core 
CPU. On the basis of MD simulations we have investigated the radial distribution function 
and the diffusion coefficient at molten and solidified states of the system. We have analysed 
the solidification kinetics of the clusters as a function of the over-cooling temperature and 
shown that the kinetics of the phase transition can be described within the framework of 
precursor formation theoretical model. Based on that theoretical model we had derived 




1/3 



(18) 
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various characteristics of the systems such as solid-liquid surface tension coefficient, rate 
of the precursors formation, phase transition temperature for the clusters of arbitrary size. 
This work highlights the computational advantages one can achieve with the use of GPUs 
and demonstrates a recipe to evaluate various thermodynamic characteristics of the finite 
system from a single set of MD simulations data. Further work can be devoted to the 
optimisation of the GPU code, and investigating the large systems consisting of mixture of 
different metals (nanoalloys). With minor modifications the code can be adopted for large- 
scale simulations of nanoindentation processes in various metal alloys and investigation of 
dislocation dynamics. 



ACKNOWLEDGMENTS 



The authors acknowledge Frankfurt Center for Scientific Computing for the possibility 
to perform complex computer simulation using CPU and GPU. A.V.Y. thanks Stiftung 
Polytechnische Gesellschaft Frankfurt am Main for financial support. 



1 A. Aguado and M. Jarrold, Annu. Rev. Phys. Chem. 62, 151 (2011). 

2 D. Turnbull, J. Appl. Phys. 21, 1022 (1950). 

3 D. Turnbull, J. Comp. Phys. 20, 411 (1952). 

4 P. Pawlow, Z. Phys. Chem. 65, 545 (1909). 

5 P. Buffat and J.-P. Borel, Phys. Rev. A 13, 2287 (1976). 

6 T. Castro, R. Reifenberger, E. Choi, and R. Anders, Phys. Rev. B 42, 8548 (1990). 

7 S. Lai, J. Guo, V. Petrova, G. Ramanath, and L. Allen, Phys. Rev. Lett. 77, 99 (1996). 

8 C. Bottani, A. Bassi, and B. Tanner, Phys. Rev. B 59, R15601 (1999). 

9 M. Schmidt, R. Kusche, W. Kronmiiller, B. von Issendorff, and H. Haberland, Phys. Rev. Lett. 
79, 99 (1997). 

M. Schmidt, R. Kusche, B. von Issendorff, and H. Haberland, Nature 393, 238 (1998). 

1 R. Kusche, T. Hippler, M. Schmidt, B. von Issendorff, and H. Haberland, Eur. Phys. J. D 9, 
1 (1999). 

22 



12 H. Haberland, T. Hippler, J. Donges, O. Kostko, M. Schmidt, B. von Issendorff, and H. Haber- 
land, Eur. Phys. J. D 9, 1 (1999). 

13 F. Calvo and F. Spiegelmann, J. Comp. Phys. 112, 2888 (2006). 

14 J. Reyes-Nava, I. Garzon, and K. Michaelian, Phys. Rev. B 67, 165401 (2003). 

15 F. Calvo and F. Spiegelmann, J. Comp. Phys. 120, 9684 (2004). 

16 K. Manninen, A. Rytkonen, and M. Manninen, Eur. Phys. J. D 29, 39 (2004). 

17 S. Chacko, D. Kanhere, and S. Blundell, Phys. Rev. B 71, 155407 (2005). 

18 A. Aguado and J. Lopez, Phys. Rev. Lett. 94, 233401 (2005). 

19 A. Aguado, J. Phys. Chem. B 109, 13043 (2005). 

20 E. Noya, J. Doye, D. Wales, and A. Aguado, Eur. Phys. J. D 43, 57 (2007). 

21 A. Shvartsburg and M. Jarrold, Phys. Rev. Lett. 85, 2530 (2000). 

22 G. Breaux, R. Benirschke, T. Sugai, B. Kinnear, and M. Jarrold, Phys. Rev. Lett. 91, 215508 
(2003). 

23 A. Lyalin, A. Hussien, A. Solovyov, and W. Greiner, Phys. Rev. B 79, 165403 (2009). 

24 F. Ding, K. Bolton, and A. Rosn, J. Vac. Sci. Technol. A 22, 1471 (2004). 

25 A. Jiang, N. Awasthi, A. Kolmogorov, W. Setyawan, A. Borjesson, K. Bolton, A. Harutyunyan, 
and S. Curtarolo, Phys. Rev. B 75, 205426 (2007). 

26 O. Petrov and I. Furo, Phys. Rev. E 73 (2006). 

27 A. V. Ruban and I. A. Abrikosov, Reports on Progress in Physics 71, 046501 (2008). 

28 A. Aguado, J. M. Lopez, J. A. Alonso, and M. J. Stott, J. Phys. Chem. B 105, 2386 (2001). 

29 C. M. Goringe, D. R. Bowler, and E. Hernandez, Reports on Progress in Physics 60, 1447 
(1999). 

30 Q. Yue, T. Cagin, W. Johnson, and W. Goddard, J. Comp. Phys. 115, 385 (2001). 

31 Y. Qi, T. Cagin, Y. Kimura, and W. A. Goddard, Phys. Rev. B 59, 3527 (1999). 

32 Y. Shibuta and T. Suzuki, Chem. Phys. Lett. 498, 323 (2010). 

33 Y. Shibuta and T. Suzuki, Chem. Phys. Lett. 502, 82 (2011). 

34 W. G. Ill, D. Brenner, S. Lyshevski, and G. I. (eds.), Handbook of Nanoscience, Engineering, 
and Technology (CRC Press, 2007). 

35 A. Harutyunyan, T. Tokune, and E. Mora, Appl. Phys. Lett. 87, 051919 (2005). 

36 O. Obolensky, I. Solov'yov, A. Solov'yov, and W. Greiner, in International Symposium "Atomic 
Cluster Collisions: structure and dynamics from the nuclear to the biological scale", Vol. 31D, 

23 



edited by A. Solov'yov (European Physical Society, 2007) p. 176. 

37 A. Harutyunyan, N. Awasthi, A. Jiang, W. Setyawan, E. Mora, T. Tokune, K. Bolton, and 
S. Curtarolo, Phys. Rev. Lett. 100, 195502 (2008). 

38 I. Solov'yov, A. Yakubovich, P. Nikolaev, I. Volkovets, and A. Solov'yov, 



J. Comp. Chem., doi 10.1002/jcc.23086 (2012)[ 



39 S. Nayak, S. Khanna, B. Rao, and P. Jena, J. Phys. Chem. A 101, 1072 (1997). 

40 M. Daw and M. Baskes, Phys. Rev. Lett. 50, 1285 (1983). 

41 M. Daw and M. Baskes, Phys. Rev. B 29, 6443 (1984). 

42 M. Finnis and J. Sinclair, Philos. Mag. 50, 45 (1984). 

43 S. Foiles, M. Daw, and M. Baskes, Phys. Rev. B 33, 7983 (1986). 

44 A. Sutton and J. Chen, Philos. Mag. Lett. 61, 139 (1990). 

45 A. Sutton, P. Godwin, and A. Horsfield, MRS Bull. 21, 42 (1996). 

46 H. Raffi-Tabar and A. Sutton, Philos. Mag. Lett. 63, 217 (1991). 

47 B. Todd and R. Lynden-Bell, Surf. Sci. 287, 191 (1993). 

48 R. Lynden-Bell, J. Phys.: Condens. Matter 7, 4603 (1995). 

49 J. Doye and D. Wales, New J. Chem. 22, 733 (1998). 

50 L. Landau and E. Lifshitz, Statistical physics, Part I (Butterworth-Heinemann, Oxford, 1980). 

51 S. Kirkpatrick, C. Gelatt, and M. Vecchi, Science 220, 671 (1983). 

52 V. Dick, I. Solov'yov, and A. Solov'yov, Phys. Rev. B 84, 115408 (2011). 

53 Y. Shibuta and T. Suzuki, Chem. Phys. Lett. 498, 323 (2010). 



24 



